Modeling of Spiking-Bursting Neural Behavior Using Two-Dimensional Map 
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A simple model that replicates the dynamics of spiking and spiking-bursting activity of real bio- 
logical neurons is proposed. The model is a two-dimensional map which contains one fast and one 
slow variable. The mechanisms behind generation of spikes, bursts of spikes, and restructuring of 
the map behavior are explained using phase portrait analysis. The dynamics of two coupled maps 
which model the behavior of two electrically coupled neurons is discussed. Synchronization regimes 
for spiking and bursting activity of these maps are studied as a function of coupling strength. It is 
demonstrated that the results of this model are in agreement with the synchronization of chaotic 
spiking-bursting behavior experimentally found in real biological neurons. 
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I. INTRODUCTION 

Understanding dynamical principles and mechanisms 
behind the control of activity, signal and information pro- 
cessing that occur in neurobiological networks is hardly 
possible without numerical studies of collective dynamics 
of large networks of neurons. These simulations need to 
take into account the complex architecture of couplings 
among individual neurons that is suggested by data from 
biological experiments. One of the complicating factors 
in understanding the simulation results is the complex- 
ity of temporal behavior of individual biological neurons. 
This complexity is due to the large number of ionic cur- 
rents involved in the nonlinear dynamics of neurons. As 
a result, realistic channel-based models proposed for a 
single neuron is usually a system of many nonlinear dif- 
ferential equations (see, for example and the re- 
view of the models in |^). The strong nonlinearity and 
high dimensionality of the phase space is a significant 
obstacle in understanding the collective behavior of such 
dynamical systems. The dynamical mechanisms respon- 
sible for restructuring collective behavior of a network 
of channel-based neuron models are difficult or impossi- 
ble to analyze because these mechanisms are well hid- 
den behind the complexity of the equations. If there is 
a chance to identify possible dynamical mechanisms be- 
hind a particular behavior of the network without the use 
of complex models, then this knowledge can be used to 
guide through numerical study of this behavior with the 
high-dimensional channel-based models. Based on these 
results one can specify the actual dynamical mechanism 
occurred in the network and understand the biological 
processes which contribute to it. 

One way to reveal the dynamical mechanisms is to 
use a simplified or phenomenological model of the neu- 
ron. However, many experiments indicate the restruc- 
turing of collective behavior utilizes a variety of dy- 
namical regimes generated by individual dynamics of 
a neuron. When neural dynamics include the regime 



chaotic spiking-bursting oscillations the system of dif- 
ferential equations describing each neuron should be at 
least a three-dimensional system (see for example 0). 
Further simplification of phenomenological models for 
complex dynamics of the neuron can be obtain using 
dynamical systems in the form of a map. Despite the 
low-dimensional phase space of this nonlinear map, it is 
able to demonstrate a large variety of complex dynamical 
regimes. 

The use of low-dimensional model maps can be use- 
ful for understanding the dynamical mechanisms if they 
mimic the dynamics of oscillations observed in real neu- 
rons, show correct restructuring of collective behavior, 
and are simple enough to study the reasons behind such 
restructuring. The development of a model map which is 
capable of describing various types of neural activity, in- 
cluding generation of tonic spikes, irregular spiking, and 
both regular and irregular bursts of spikes, is the goal of 
this paper. 

It is known that constructing a low-dimensional system 
of differential equation which is capable of generating fast 
spikes bursts excited on top of the slow oscillations, one 
needs to consider a system which has both slow and fast 
dynamics (see for example Using the same ap- 

proach one can construct a two-dimensional map, which 
can be written in the form 



Xn+l = f{Xn,yn + Pn) , 
Vn+l = Vn - t^{Xn + 1) + MC^n 



(la) 

(lb) 



where Xn is the fast and j/„ is the slow dynamical vari- 
able. Slow time evolution of y„ is due to small values of 
the parameter /i = 0.001. Terms /3„ and cr„ describe ex- 
ternal influences applied to the map. These terms model 
the dynamics of the neuron under the action of the exter- 
nal DC bias current Idc and synaptic inputs X^*'". The 
term (t„ can also be used as the control parameter to 
select the regime of individual behavior. 

A model in the form of two-dimensional map, whose 
form is similar to (1) was used in the study of dynami- 
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cal mechanisms behind the emergence and regularization 
of chaotic bursts in a group of synchronously bursting 

The effect of 



cells coupled through a mean field |12 



anti-phase regularization was also modeled before with 
one-dimensional maps . In both these models the os- 
cillations during the burst were described by a chaotic 
trajectory. Map (f) improves these models by adding a 
feature that enables one to mimic the dynamics of in- 
dividual spikes within the burst. This is achieved us- 
ing a modification of the shape of the nonlinear function 
f{x, y) which is now a discontinuous function of the form 




x < 

< x < a - 

X > a + y 



(2) 



where a is a control parameter of the map. The depen- 
dence of /(x, v) on X computed for a fixed value of y is 
shown in Fig.|l|. In this plot the values of a and y are 
set to illustrate the possibility of coexistence of limit cy- 
cle, Pfc, corresponding to spiking oscillation in (|la|), and 
fixed points Xps and Note that when y increases or 
decreases the graph of f{x, y) moves up or down, respec- 
tively, except for the third interval a; > a + y, where the 
values of /(x, y) always remain equal to -1. 




FIG. 1. The shape of the nonlinear function f{x,y), plot- 
ted for a = 6.0 and y — —3.93, is shown with a solid line. 
The dashed line illustrates a super-stable cycle Pk of the fast 
map (|3a|), where the value of y is fixed. The stable and un- 
stable fixed points of the map are indicated by Xps and Xpu, 
respectively. 

The paper is organized as follows: Section |l| considers 
the features of the fast and slow dynamics which explain 
the formation of various types of behavior in the isolated 
map. It also discusses the bifurcations responsible for the 
qualitative change of activity. The dynamics of response 
generated by the map after it is excited or inhibited by 
an external pulse is discussed in Section III . The goal of 



this section is to illustrate dynamical mechanisms that 
control the response, and to show how the relation be- 
tween the two inputs of the map influence the properties 



of the response. Section IV presents the results of a syn- 



chronization study in two coupled maps. 



II. INDIVIDUAL DYNAMICS OF THE MAP 

Consider the regimes of oscillations produced by the 
individual dynamics of the map. In this case the inputs 
(3n = (3 and an = o are constants. Note that if /? is a 
constant, then it can be omitted from the equations us- 
ing the change of variable yn+ P ^ ■ Therefore, the 
individual dynamics of the map depends only on the con- 
trol parameters a and cr, and the map can be rewritten 
in the form 



yn+1 =yn- t^{Xn + I) + tl(7 



(3a) 
(3b) 



Typical regimes of temporal behavior of the map are 
shown in Figs || and |^. When the value of a is less then 
4.0 then, depending on the value of parameter cr, the map 
generates spikes or stays in a steady state (see Fig^). The 
frequency of the spikes increases as the value of parameter 
(T is increased (see Fig ^) . 
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FIG. 2. Waveforms of spiking behavior generated by map 
(3) with a = 4. Panel (a) shows the transition to the regime 
of silence for cr = —0.01. The regimes of continuous tonic 
spiking are computed for o — 0.01 (b) and a — 0.1 (c). 

For a > 4 the map dynamics are capable of producing 
bursts of spikes. The spiking- bursting regimes are found 
in the intermediate region of parameter a between the 
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regimes of continuous tonic spiking and steady state (si- 
lence). The spiking-bursting regimes include both peri- 
odic and chaotic bursting. A few typical bursting-spiking 
regimes computed for different values of a are presented 
in Fig I 

Due to the two different time scales involved in the 
dynamics of the model, the mechanisms behind the re- 
structuring of the dynamical behavior can be understood 
through the analysis of the fast and slow dynamics sepa- 
rately. In this case, time evolution of the fast variable, x, 
is studied with the one-dimensional map ( pa| ) where the 
slow variable, ?/, is treated as a control parameter whose 
value drifts slowly in accordance with equation (^). It 
follows from (3b) that the value of y remains unchanged 
only '\i X — Xs given by 



X. = -1 



(4) 



\i X < Xs, then the value of y slowly increases. \i x > Xs-, 
then y decreases. 
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FIG. 3. Typical waveforms of the spiking-bursting behavior 
generated by the map (3) for the following parameter values: 
a = 4.5, G = 0.14 (a); a = 6.0, a = -0.1 (b); and a = 6.0, 
a = 0.386 (c). 



From ( 3a ) one can find the equation for the coordinate 
of fixed points Xp of the fast map. This equation is of 
form 



y 



1 X71 



(5) 



where Xp < 0. Equation (||) defines the branches of slow 
motion in the two-dimensional phase space {xn,yn) (see 
Fig H). The stable branch Sps{y) exists for Xp < I — %foL 
and the unstable branch Spuiy) exists within 1 — ^/a < 



< 



Considering the fast and slow dynamics together, one 
can see that, if Xs is in the stable branch Sps{y), then 
the map (3) has a stable fixed point. The stable fixed 
point corresponds to the regime of silence in the neu- 
ral dynamics. The oscillations in the map dynamics will 
appear when Xg > I — y/ct- This is the threshold of exci- 
tation, which corresponds to the bifurcation values of a 
given by 



(6) 



To understand the dynamics of excited neurons within 
this model one needs to consider branches that corre- 
spond to the spiking regime. To evaluate the location 
of the spiking branches, consider the mean value of a;„ 
computed for the periodic trajectory of the fast map (|3a| ) 
as a function of y. Here y is treated as a parameter. Use 
of such an approximation is quite typical for an analysis 
involving fast and slow dynamics. It works well for small 
values of parameter ^. 

It follows from the shape of f{x,y) that for any value 
of y, map ( ^a|) generates no more than one periodic tra- 
jectory Pfc, where k is the period (i.e. a;„ — Xn+k)- The 
cycles Pk always contain the point x^ — —1 and the 
index k increases stepwise {k ^ k + 1) as y decreases. 
Since the trajectory always has a point in the flat in- 
terval of f{x,y), all these cycles are super-stable except 
for bifurcation values of y for which the trajectory con- 
tains the point Xn = (see Fig |l]) . The location of the 
spiking branch, S spikes, of "slow" motion in the phase 
plane {xn,yn) can be estimated as the mean value of x 
computed for the period of cycle Pk, 
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k 

k ^ 

m— 1 



(7) 



where k is the period of Pk, and f'^"^\x,y) is the m-th 
iterate of (^) , started at point x and computed for fixed 
y. The spiking branch of "slow" dynamics evaluated with 
(0) is shown in Fi^. One can see that this branch has 
many discontinuities caused by the bifurcations of the 
super-stable cycles Pk- 

To complete the picture of fast and slow dynamics of 
the model for a > 4, one needs to consider the fast map 
bifurcation associated with the formation of homoclinic 
orbit hpu originating from the unstable fixed point, Xpu- 
This homoclinic orbit occurs when the coordinate of Xpu 
become equal to -1. It can be easily shown that such a 
situation can take place only if a > 4. The homoclinic 
orbit forms at the value of y where the unstable branch 
Spuiy) crosses the line x = —1, see Fig.^. When the 
map is firing spikes and the value of y gets to the bi- 
furcation point, the cycle Pk merges into the homoclinic 
orbit, disappears, and then the trajectory of the map 
jumps to the stable fixed point Xpg. 

Typical phase portraits of the model, obtained un- 
der the assumptions made above, are presented in Figj^. 
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Fig.^ shows the typical behavior for 2 < a < 4. Here, 
only two regimes are generated. The first regime is the 
state of silence, when the operating point (OP), given by 
the intersection of = —1 + a and one of the branches, 
is on the branch Sps (y) . The second regime is the regime 
of tonic spiking, when the operating point is selected on 
the spiking branch S spikes ■ 
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FIG. 4. Stable {Sps{y), Sspikes) and unstable {Spu{y)) 
branches of slow dynamics of (3) plotted on the plane of phase 
variables {y,x). The cases of a = 4 and a = 6 are presented 
in (a) and (b), respectively. Numbers on the branch Sspikes 
stand for the value of k of the cycles Pk- The operating points 
OP are selected to illustrate the regime of silence in (a) and 
the regime of spiking-bursting oscillations in (b). Arrows in- 
dicate the direction of slow evolution along the branches, and 
switching in the case of (b). 

When a > 4 the picture changes qualitatively (see 
Fig ^). Now, stable branches Sps{y) and Sspikes are 
separated by the unstable branch Spu{y)- If the operat- 
ing point is selected on Spu{y), then the phase of silence, 
corresponding to slow motion along Sps{y), and spiking, 
when system moves along Sspikes, alternate forming the 
regime of spiking-bursting oscillations. The beginning of 
a burst of spikes corresponds to the bifurcation state of 
the fast map where fixed points Xps and Xpu merge to- 



gether and disappear. Before this bifurcation the system 
is in Xps and, therefore, y increases. The termination of 
the burst is due to the bifurcation of the fast map as- 
sociated with the formation of the homoclinic orbit hp^. 
Here, the limit cycle of the spiking mode merges into 
the homoclinic orbit and disappears. After that the fast 
subsystem flips to the stable fixed point Xps- Then the 
process repeats (see arrows in Fig ||d). 

It is clear from Figji^ that when operating point is 
set on Sps{y) the model will be in the regime of silence. 
When the operating point is set on the branch Sspikes 
the model produces tonic spiking, unless the point is set 
close to the formation of homoclinic orbit hpu- One can 
see that, at the vicinity of this bifurcation, the branch 
Sspikes becomes densely folded. As a result, the behav- 
ior of y, which is governed by the mean value of Xn, 
can become extremely sensitive to small perturbations 
and even lead to instability caused by high-gain feed- 
back. This is one of the reasons for the irregular, chaotic 
spiking-bursting behavior which occurs in the map when 
the operating point is set close to the area of the branch 
Sspikes where this branch is densely folded. The detailed 
and rigorous analysis of chaotic dynamics cannot be done 
within the approximations made above and require more 
precise computation of Sspikes which is beyond the scope 
of this paper. 
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FIG. 5. Bifurcation diagram on the parameter plane ((t, a). 

The results of the analysis presented above are sum- 
marized in the sketch of the bifurcation diagram plotted 
on the parameter plain (cr, a) (see Fig ||) . The bifurca- 
tion curve ath corresponds to the excitation threshold (|^) 
where the fixed point of the 2-d map becomes unstable 
and the map starts generating spikes. Curve Lts shows 
the approximate location of the border between spiking 
and spiking-bursting regimes, obtained in the numerical 
simulations. Note that separation of spiking and burst- 
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ing regimes is not always obvious, especially in the regime 
of chaotic spiking. The regime of spiking-bursting oscil- 
lations takes place within the upper triangle formed by 
curves ath and Lta- The regimes of chaotic spiking or 
spiking-bursting behavior are found in the relatively nar- 
row region of the parameters located around Lts- This 
region contains complex structure of bifurcations, associ- 
ated with multi-stable regimes, and is not shown in Fig|^. 

The bifurcation diagram shows the role of control pa- 
rameters in the selection of dynamical features of the 
considered neuron model. Both parameters can be used 
to mimic a particular type of neural behavior. Parameter 
(J can be used to model the external DC current injec- 
tion that depolarizes or hyperpolarizes the neuron. In 
this case a can be written as 



a — 



'DC, 



(8) 



where is the parameter that selects the dynamics of 
isolated neuron, and Idc is the parameter that models 
the DC current injected into the cell. The changes in be- 
havior caused by Idc are similar to the action of the pa- 
rameter / in the well known Hindmarsh-Rose model |Q. 

If the individual dynamics of the modelled neuron are 
capable of generating only regimes of silence and tonic 
spiking, and do not support the regime of bursts of spikes, 
then the value of a should be set below 4. In this case, 
no matter what the level of Idc injected is, the system 
will not show bursts of spikes (see Fig ^). 



III. MODELING OF RESPONSE TO THE 
INJECTION OF CURRENT 

To study the dynamical regimes of neuronal behavior 
in experiments, biologists change the type of neural ac- 
tivity using the injection of electrical current into the cell 
though an electrode. It was shown above that the injec- 
tion of DC current can be modeled in the map (1) using 
parameter cr„ = cr (see equation (^)). In this case, since 
the external influence does not vary in time, the role of 
parameter /?„ = /3 is not important, because the behav- 
ior of the map after the transient is independent of the 
value of (3. However, when the injected current changes 
in time, it may be useful to consider the dynamics of 
parameter /3„ in order to provide more realistic modeled 
behavior during the transient. Taking this into account, 
the input of the model can be considered in the form: 



Pn = P'ln 



(9) 



where /„ is injected current, and coefficients P'^ and a'^ 
are selected to achieve the desired properties of response 
behavior. 

This Section briefly illustrates how the relation be- 
tween these coefficients effects the dynamics of response 
to the pulse of /„. To be specific, consider the map in the 



regime of tonic spiking with a = 5.0, a = 0.33, /3 = 0. To 
study the response behavior, positive and negative pulses 
of amplitude 0.8 and duration of 100 iterations were ap- 
plied to the continuously spiking map. In the simula- 
tions presented below the coefficient a'^ was selected to 
be equal to one. 
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FIG. 6. Response of the model with a'^ = 1 and = to 
a positive pulse of /„. The parameters of the map are selected 
in the regime of tonic spiking (a — 5.0, a = 0.33, /3 = 0). 
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FIG. 7. Response of the model with cr' = 1.0 and (3" — 
to a negative pulse of /„. The parameters of the map are 
selected in the regime of tonic spiking [a — 5.0, a — 0.33, 
13 = 0). 

Figures |^ and ^ show the response to the positive and 
negative pulse, respectively, computed with P*^ = 0. In 
this case, during the action of a positive pulse, the value 
of Un increases monotonically because of the increased 
value of (T„ (see ([lb[)). The increase of ?/„ pushes the 
fast map up, see Fig ^ This leads to an increase in 
the frequency of spiking. After the action of the pulse 
ends the value of j/„ monotonically decreases back to the 
original state (see Fig. ||). 

When a negative pulse is applied, it pushes the operat- 
ing point down and, if the amplitude of the pulse is suf- 



5 



ficiently large, shuts off the regime of spiking (see Fig.^. 
This happens because the trajectory of the fast map gets 
to the perturbed operating point which is now on the 
stable branch Sps{y), (see Fig ^d). After the pulse is 
over the spiking does not re-appear immediately because 
the system spends some time drifting along the stable 
branch of slow motions Sps (y) , during which the variable 
y overshoots its original value for the spiking regime. As 
a result, after the system switches to the spiking branch 
S spikes 1 Vn monotonically drifts down to the unperturbed 
operating point. The dynamics of the slow evolution are 
clearly seen in the lower panel of FigJ^. 
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FIG. 8. Response of the model with a" = 1.0 and (3" = 1.0 
to a positive pulse of /„. The parameters of the map are 
selected in the regime of tonic spiking (a = 5.0, a — 0.33, 
13 = 0). 
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FIG. 9. Response of the model with a" = 1.0 and 13" = 1.0 
to a negative pulse of I„. The parameters of the map are 
selected in the regime of tonic spiking {a = 5.0, a — 0.33, 
13 = 0). 

Figures |^ and ^ illustrate how response to the pulse 
of /„ changes when coefficient /3„ is not equal to zero. 
To be brief and specific consider the case of P'^ = 1.0 
and cr'^ = 1.0. Comparing FigJ3 with Fig.pl one can see 



that the response dynamics to the positive pulse changes 
qualitatively. Indeed, when the pulse current is applied, 
then, acting through /3„, it immediately forces the fast 
map to shift up. As a result, the rate of spiking and 
mean value of Xn increases sharply. Variable y„ reacts to 
this change and decreases its value to compensate for the 
sudden change of the mean value of x„. When the action 
of the pulse is over, the value of /3n returns to its original 
value and, due to the updated levels of yn, the fast map 
overshoots its original state. As a result, the trajectory 
of the fast map reaches the stable fixed point. To return 
to the original regime of spiking the system has to go all 
the way along the branches of slow motion Spgijj) and 
S spikes back to the original operating point (see Fig. ^). 
This type of response is observed in real neurobiological 
experiments, see for example [ p^ . 

Using the same analysis one can understand the new 
effects in the response to a negative pulse caused by the 
action of /3„. These new effects can be clearly seen com- 
paring Fig.^ with Fig. |^. 

The presented results illustrate how 2-d map (1) along 
with equation (^ can model a large variety of transient 
neural behavior induced by injected current. Due to the 
simplicity of the model one can clearly see the nonlin- 
ear mechanisms behind the response behavior and apply 
them to select the desired balance between cr^ and (3'^ to 
model a particular type of response. 

IV. REGIMES OF SYNCHRONIZATION IN TWO 
COUPLED MAPS 

This section presents the results of studies of syn- 
chronization regimes in coupled chaotically bursting 2-d 
maps. The goal of this study is to reproduce the main 
regimes of synchronous behavior found in a real neuro- 
biological experiment ||l^ . The experiment was carried 
out on two electrically coupled neurons (the pyloric dila- 
tors, PD) from the pyloric CFG of the lobster stom- 
atogastric ganglion [p^ . The regimes found in the ex- 
periment were also reproduced in numerical simulations 
using Hindmarsh-Rose model ||l8|,|9|, one-dimensional 
map model jl^ , and in experiments with electronic neu- 
rons (l|]. 

The equations used in the numerical simulations of the 
coupled maps are of form 

-^z,n-fl — f{-^i,nj yi.n ^" l^i.ri)^ 

yi,n+l = yi,n - K^i.n + 1) + M^^i + M^'i.n, (10) 

where index i specifies the cell, and ai is the parameter 
that defines the dynamics of the uncoupled cell. The cou- 
pling between the cells is provided by the current flowing 
from one cell to the other. This coupling is modeled by 

Pi,7i — 9jiP {"^j^n •^i,n) 
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where i j, and gji is the parameter characterizing the 
strength of the couphng. The coefficients f]'^ , set the 
balance between the coupHngs for the fast and slow pro- 
cesses in the cells, respectively. In the numerical simu- 
lations the values of the coefficients are set to be equal: 
(3^ = 1.0 and a'^ ~ 1.0. The other parameters of the 
coupled maps (|lO|) that remain unchanged in the simu- 
lations have the following values: n = 0.001, ai — 4.9, 
a2 = 5.0. The coupling between the maps is symmetri- 
cal, = gij = g. 
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FIG. 10. Waveforms generated by the two coupled 2-d 
maps with cti = 0.240 and (T2 = 0.245. Waveforms of 
a;i,„ and X2,n are shown by solid and dashed lines, respec- 
tively. Three different regimes of spiking-bursting behavior 
are shown. Individual chaotic spiking-bursting oscillations 
of uncoupled maps, g = 0.0, - (a). Regime of synchronized 
chaotic bursts, computed with g = 0.043, - (b). Regime of 
antiphase synchronization, computed with g = —0.029, -(c). 

First, consider the main regimes of synchronization be- 
tween the maps generating irregular, chaotic bursts of 
spikes. To set the individual dynamics of the maps to 
this regime, the parameters ai, that take into account 
the DC bias current injected into the neurons (see |jl^ for 
details), are tuned to the following values: ai = 0.240, 
a2 — 0.245. When these maps are uncoupled {g — 0) 
they produce chaotic spiking-bursting oscillations shown 



in Fig 10 a,. When the coupling becomes sufficiently large 
the slow components of the bursts synchronize while 
spikes within the bursts remain asynchronous see the 
waveforms in Fig. obtained with g — 0.043. This 
regime of synchronization is typical for naturally coupled 
PD neurons (see Figure 2a in |p^ ). 

Introduction of negative coupling, g < 0, also leads to 
synchronization of bursts, but in this regime of synchro- 
nization the systems burst in antiphase. Typical wave- 
forms produced in this regime are presented in Fig. |lC~ 



This regime of synchronization of chaotic bursts is also 
observed in the experimental study of coupled PD neu- 
rons (see Figure 2c in [^). It is important to emphasize 
that, both in this simulation and in the experiment, the 
regime of antiphase synchronization characterized by the 
onset of regular bursts. 

The simplicity of this model enables one to understand 
the possible cause for the onset of regular bursting in the 
antiphase synchronization. It is shown in Section that 
chaotic dynamics of bursts occurs when the operating 
point of the system appears close to the leftmost area 
of the spiking branch Sspikes- In this region, Sapikes is 
densely folded and, if system slows down in this area, the 
timing for the end of the burst becomes very sensitive to 
infinitesimal perturbations. 
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FIG. 11. Attractors computed from the waveforms of the 
first cell operating in the regime of uncoupled oscillations - (a) 
and in the regime of antiphase synchronization - (b). Wave- 
forms for these regimes are shown in Fig. |lo| a and Fig. px|c, 
respectively. The points of the attractors are connected by 
lines to clearly show their sequence in time. 



This fact is illustrated in Fig. |ll|a. This figure shows 
the attractor computed from the waveforms a;i,„ and 
of the first cell when the cells are uncoupled (this regime 
is shown in Fig. p^). To plot the attractor, the wave- 
forms of xi^n and were filtered by a fourth order 
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low-pass filter with cutoff frequency 0.6. As a result, 
the attractor does not contain sharp spikes and the ap- 
proximate location S spikes, which is close to the center 
of filtered spiking oscillations, is easy to see. The level 
of coordinate x/, corresponding to the operating point, 
is equal to Xs — — 1 + cti —0.76. One can see from 
Fig. [ll| a that when the center of oscillations gets close 
to that level the trajectory becomes rather complex, as 
shown by the dense set of trajectories in the left part of 
the attractor. 

The numerical simulations show that this effect of reg- 
ularization takes place even when a'^ = 0. Therefore, the 
vertical shift of operating point, Xs, caused by the slow 
coupling current is not very important for this effect. The 
coupling term /3i.„ directly influences the fast dynamics 
by shifting the graph of the 1-d map (see Fig. p up or 
down depending on the sign of /3i,„. In the regime of 
antiphase synchronization f3i^n is positive, when the i-th 
cell is spiking, and negative, when it is silent. Therefore, 
from the viewpoint of fast dynamics, the coupling forces 
the cells to stay on the current branch of slow motion. 
One can understand these dynamics from the graph of 
/(x, y) and (p^. This effect is also seen as the formation 
of extended shape of attractor plotted for the regime of 
antiphase synchronization (see Fig. |ll|b). Note that axes 
in Fig. [ri| a and Fig. |ll|b have different scales. 




100 200 300 400 500 

n 



FIG. 12. Tonic spiking waveforms generated with 
CTi = 0.653 and (T2 = 0.714. Waveforms of xi^n and X2,n 
are shown by solid and dashed lines, respectively. Spiking in 
the uncoupled maps, g — 0.0, - (a). Regime of synchronized 
spikes, computed with g = 0.008, - (b) 

In the regime of chaotic bursts the level of Xg is close to 
the stable branch of spiking Sspikes, and, therefore, the 
slow evolution along the branch of silence Spu is faster 
then along Sspikes ■ This means that duration of the phase 
of silence in the cell is shorter than duration of the burst 
of spikes. When the cell switches from silence to spiking 
it sharply changes the value and the sign of This 
change pushes Sspikes of the spiking cell down to its origi- 
nal location, and, as a result, quickly drags the trajectory 



through the area of complex behavior toward the branch 
Spu- Due to this fast change the duration of bursting 
is determined only by the dynamics of the silent phase, 
which is regular. 

It was observed in the experiment that, in the regime 
of tonic spiking, spikes in PD neurons synchronize at low 
levels of coupling (i.e. without added artificial coupling, 
see ||l^ for details). This property of synchronization in 
the regime of tonic spiking is also typical for the dynam- 
ics of the maps considered here (see Fig.^2|). In these 
numerical simulations, the uncoupled maps were tuned 
to generate spikes with slightly different rates. One can 
clearly see the beats between the waveforms plotted in 
Fig.|l^a. When small coupling, g = 0.008, is introduced 
the beats disappear as soon as the spikes get synchro- 
nized (see Fig. [l^). 

Although the synchronization between continuously 
spiking maps is easily achieved, it is important to un- 
derstand that the dynamics of such synchronization in 
the maps is a much more complex process than it is in 
the models with continuous time. Due to the discrete 
time, the periodic spikes can lock with different frequency 
ratios. This results in the existence of a complex multi- 
stable structure of synchronization regimes. This struc- 
ture becomes more noticeable in the dynamics of syn- 
chronization as the number of iterations in the period of 
spikes decreases. 

V. DISCUSSION AND OUTLOOK 

The simple phenomenological model of complex dy- 
namics of spiking-bursting neural activity is proposed. 
The model is given by two-dimensional map (1). The 
first equation (|la|) of the map describes the fast dynam- 
ics. Its isolated dynamics is capable of generating sta- 
ble limit cycles, which mimics the spiking activity of a 
neuron, and a stable fixed point, which corresponds to 
phase of silence. This 1-d map has a region of param- 
eters were both these stable regimes coexist. Existence 
of such multistable regimes is the reason for generation 
of bursts, when the operating point, defined by the dy- 
namics of the second equation (p^), is properly set (see 
Fig Id). 

The shape of the nonlinear function f{x,y) used in 
(|la| ) is selected in the form (^ based on the following 
considerations: (i) The fast map should generate limit 
cycles whose waveforms mimics those of the spikes, (ii) 
Each spike generated by the map always has a single 
iteration residing on the most right interval x > a + 
y. Therefore, the moment of time corresponding to the 
appearance of a trajectory on this interval can be used to 
define the time of a spike. This feature is important for 
modeling the dynamics of chemical synapses in a group 
of coupled neurons, (iii) The analytic expression of the 
nonlinear function in interval a; < should be simple 
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enough to allow rigorous analysis of bifurcations of fixed 
points of the fast map. (iv) Use of the fixed level, —1, 
in the rightmost interval of the function simplifies the 
analysis of dynamics at the end of the bursts. The end of 
a burst is associated with the formation of a homoclinic 
orbit, which corresponds to the case when the trajectory 
from this interval maps into the unstable fixed point (see 
Section H). 

It is clear that the shape of function (||) can be modified 
to take into account other dynamical features that need 
to be modeled. Due to the low-dimensionality of the 
model, the dynamical mechanisms behind its behavior 
are easy to understand using phase plane analysis. This 
allows one to see how introduced modifications influence 
the dynamics. 

Some modification of the fast map is required to en- 
hance the region of parameters where the model can still 
be used to mimic neural dynamics properties. For ex- 
ample, from (^) and Fig. |l| one can see that, due to the 
dynamics of coupling terms, the trajectory of the fast 
system can stay in the middle interval of f{x, y) for sev- 
eral iterations. This can happen when the external in- 
fluence monotonically pushes the function up while the 
trajectory is located in the middle interval. In this case 
the trajectory will map to the middle interval again and 
again, increasing the duration of a spike. This artifact 
can be removed using one of the following modifications. 
Introduction of a sufficient gap between the right end of 
this interval and the diagonal (see Fig. ^ will help the 
map to terminate the spike despite the monotonic ele- 
vation of the function. Alternatively one can introduce 
an additional condition to the fast map ( p^ ) that forces 
the map always to iterate its trajectory from the mid- 
dle interval to the rightmost one, despite the dynamics 
of y (see (||)). This can be achieved using the function 
f{xmy) of the following form 

{a/{l-x„) + y, Xn<0 
a + y, 0<Xn<a + y 

— 1, a;„ > a + ?/ or > 

An important feature of the model discussed here is 
that one can use two inputs, /3„ and cr„, to achieve the 
desired response dynamics. Although these inputs are 
not directly related to dynamics of specific ionic currents, 
they can be used to capture the collective dynamics of 
these currents. Selecting a proper balance between these 
inputs, one can model a large variety of the responses 
that are seen in different neurons. Again, the simplic- 
ity of the model helps one to understand the dynamical 
properties of each input and set the proper balance be- 
tween them. 
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